Compare urban traits and phylogency

source('../env.R')

Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
── Attaching core tidyverse packages ──────────────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.0.2     ── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errorsRegistered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      

Attaching package: ‘dbplyr’

The following objects are masked from ‘package:dplyr’:

    ident, sql

Loading required package: ape

Attaching package: ‘ape’

The following object is masked from ‘package:dplyr’:

    where

Loading required package: maps

Attaching package: ‘maps’

The following object is masked from ‘package:purrr’:

    map

Please cite the eBird Status and Trends data using: 
  Fink, D., T. Auer, A. Johnston, M. Strimas-Mackey, S. Ligocki, O. Robinson, 
  W. Hochachka, L. Jaromczyk, C. Crowley, K. Dunham, A. Stillman, I. Davies, 
  A. Rodewald, V. Ruiz-Gutierrez, C. Wood. 2023.
  eBird Status and Trends, Data Version: 2022; Released: 2023. Cornell Lab of
  Ornithology, Ithaca, New York. https://doi.org/10.2173/ebirdst.2022

This version of the package provides access to the 2022 version of the eBird
Status and Trends Data Products. Access to the 2022 data will be provided 
until November 2024 when it will be replaced by the 2023 data. At that 
point, you will be required to update this R package and transition to using 
the new data.
terra 1.8.5

Attaching package: ‘terra’

The following object is masked from ‘package:phytools’:

    rescale

The following objects are masked from ‘package:ape’:

    rotate, trans, zoom

The following object is masked from ‘package:tidyr’:

    extract

Loading required package: permute
Loading required package: lattice
This is vegan 2.6-8

Attaching package: ‘vegan’

The following object is masked from ‘package:phytools’:

    scores

Loading required package: mvtnorm
Loading required package: survival
Loading required package: TH.data
Loading required package: MASS

Attaching package: ‘MASS’

The following object is masked from ‘package:terra’:

    area

The following object is masked from ‘package:dplyr’:

    select


Attaching package: ‘TH.data’

The following object is masked from ‘package:MASS’:

    geyser

Registered S3 method overwritten by 'broom':
  method        from 
  nobs.multinom MuMIn

Attaching package: ‘ggpubr’

The following object is masked from ‘package:terra’:

    rotate

The following object is masked from ‘package:ape’:

    rotate

Loading required package: nlme

Attaching package: ‘nlme’

The following object is masked from ‘package:dplyr’:

    collapse

Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2

Attaching package: ‘GGally’

The following object is masked from ‘package:terra’:

    wrap


Attaching package: ‘foreach’

The following objects are masked from ‘package:purrr’:

    accumulate, when


Attaching package: ‘scales’

The following object is masked from ‘package:terra’:

    rescale

The following object is masked from ‘package:phytools’:

    rescale

The following object is masked from ‘package:purrr’:

    discard

The following object is masked from ‘package:readr’:

    col_factor


Attaching package: ‘cowplot’

The following object is masked from ‘package:ggpubr’:

    get_legend

The following object is masked from ‘package:lubridate’:

    stamp

── Attaching packages ─────────────────────────────────────────────────────────────────────────────────────────────────── treedataverse 0.0.1 ──
✔ tidytree    0.4.6      ✔ ggtree      3.14.0
✔ treeio      1.30.0     ✔ ggtreeExtra 1.16.0
── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────── treedataverse_conflicts() ──
✖ foreach::accumulate()  masks purrr::accumulate()
✖ ggtree::collapse()     masks nlme::collapse(), dplyr::collapse()
✖ scales::discard()      masks purrr::discard()
✖ treeio::drop.tip()     masks tidytree::drop.tip(), ape::drop.tip()
✖ ggtree::expand()       masks tidyr::expand()
✖ tidytree::filter()     masks dplyr::filter(), stats::filter()
✖ ggtree::flip()         masks terra::flip()
✖ .GlobalEnv::geom_map() masks ggplot2::geom_map()
✖ dbplyr::ident()        masks dplyr::ident()
✖ ggtree::inset()        masks terra::inset()
✖ tidytree::keep.tip()   masks ape::keep.tip()
✖ dplyr::lag()           masks stats::lag()
✖ maps::map()            masks purrr::map()
✖ treeio::mask()         masks terra::mask()
✖ treeio::read.newick()  masks phytools::read.newick()
✖ ggtree::rotate()       masks ggpubr::rotate(), terra::rotate(), ape::rotate()
✖ tidytree::select()     masks MASS::select(), dplyr::select()
✖ dbplyr::sql()          masks dplyr::sql()
✖ terra::trans()         masks ape::trans()
✖ foreach::when()        masks purrr::when()
✖ ape::where()           masks dplyr::where()
✖ terra::zoom()          masks ape::zoom()
Loading required package: usethis
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: ‘devtools’

The following object is masked from ‘package:permute’:

    check

Using GitHub PAT from the git credential store.
Skipping install of 'clootl' from a github remote, the SHA1 (30ac57ee) has not changed since last install.
  Use `force = TRUE` to force installation
To cite package ‘clootl’ in publications use:

  Miller E, McTavish E, Sanchez-Reyes L (2025). “clootl: Fetch and Explore the Cornell Lab of Ornithology Open Tree of Life
  Avian Phylogeny.” <https://github.com/eliotmiller/clootl>.

  McTavish E, Gerbracht J, Holder M, Iliff M, Lepage D, Rasmussen P, Redelings B, Sanchez-Reyes L, Miller E (2025). “A complete
  and dynamic tree of birds.” _Proceedings of the National Academy of Sciences_.

To see these entries in BibTeX format, use 'format(<citation>, bibtex=TRUE)', or 'toBibtex(.)'.

The current version of the Aves tree is v1.4.
        Please specify the tree and taxonomy version used when citing this R package.
        When possible, cite all the original studies supporting your tree:
        These citations are acessible using getCitations(your_tree)
Loading required package: spData
To access larger datasets in this package, install the spDataLarge package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`
Registered S3 method overwritten by 'spdep':
  method   from
  plot.mst ape 

Attaching package: ‘gridExtra’

The following object is masked from ‘package:dplyr’:

    combine


Attaching package: ‘grid’

The following object is masked from ‘package:terra’:

    depth

Spherical geometry (s2) switched off

Create dataset

Load community data

communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))

Rows: 2428 Columns: 7── Column specification ────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): city_name, ebird_species_name, seasonal, presence, origin
dbl (2): city_id, relative_abundance_proxy
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

communities

Load trait data

traits = read_csv(filename(TAXONOMY_OUTPUT_DIR, 'traits_ebird.csv'))

Rows: 332 Columns: 10── Column specification ────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): ebird_species_name, habitat, trophic_level, trophic_niche, primary_lifestyle
dbl (5): beak_width, hwi, mass, habitat_density, migration
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(traits)

species_urban_response = communities %>% group_by(ebird_species_name) %>% summarise(mean_relative_abundance_proxy = mean(relative_abundance_proxy)) %>% left_join(traits)

Joining with `by = join_by(ebird_species_name)`

species_urban_response$is_urban = species_urban_response$mean_relative_abundance_proxy > 0
species_urban_response$ebird_species_name = str_replace(species_urban_response$ebird_species_name,' ', '_')
species_urban_response

Most urban species

species_urban_response %>% filter(mean_relative_abundance_proxy > 50)

species_urban_response %>% filter(mean_relative_abundance_proxy > 40)

species_urban_response %>% filter(mean_relative_abundance_proxy > 30)

Difference amongst urban and non-urban

species_urban_response %>% group_by(is_urban) %>% summarise(
  median_beak_width = median(beak_width),
  iqr_beak_width = IQR(beak_width),
  
  median_hwi = median(hwi),
  iqr_hwi = IQR(hwi),
  
  median_mass = median(mass),
  iqr_mass = IQR(mass),
  
  n = n()
)

Compare traits between urban and non-urban species

Gape Width

ggplot(species_urban_response, aes(x = is_urban, y = beak_width)) + geom_boxplot()

wilcox.test(beak_width ~ is_urban, species_urban_response)

    Wilcoxon rank sum test with continuity correction

data:  beak_width by is_urban
W = 1525, p-value = 0.3893
alternative hypothesis: true location shift is not equal to 0

There is not a significant response for the gape width between urban species and non urban species.

HWI

ggplot(species_urban_response, aes(x = is_urban, y = hwi)) + geom_boxplot()

wilcox.test(hwi ~ is_urban, species_urban_response)

    Wilcoxon rank sum test with continuity correction

data:  hwi by is_urban
W = 1213, p-value = 0.01492
alternative hypothesis: true location shift is not equal to 0

There is a significant response between urban and non-urban species

Mass

ggplot(species_urban_response, aes(x = is_urban, y = mass)) + geom_boxplot()

wilcox.test(mass ~ is_urban, species_urban_response)

    Wilcoxon rank sum test with continuity correction

data:  mass by is_urban
W = 1825, p-value = 0.5168
alternative hypothesis: true location shift is not equal to 0

There is not a significant response between species.

Phylogeny

Phylogentic response to relative abundance?

tree = read.tree(filename(TAXONOMY_OUTPUT_DIR, 'phylogeny.tre'))
tree_pruned = ladderize(drop.tip(tree, setdiff(tree$tip.label, species_urban_response$ebird_species_name)))
ggtree(tree_pruned, layout="fan")

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

columbidae_response <- species_urban_response$mean_relative_abundance_proxy
names(columbidae_response) <- species_urban_response$ebird_species_name

Is the a phylogenetic signal for our relative abundance proxy

col_ut_phylo_signal <- phylosig(tree_pruned, columbidae_response, method="lambda", test=TRUE)
col_ut_phylo_signal

Phylogenetic signal lambda : 0.46832 
logL(lambda) : -576.839 
LR(lambda=0) : 7.45324 
P-value (based on LR test) : 0.00633224 

Plot tree

species_urban_response$ebird_species_name_clean = gsub("_", " ", species_urban_response$ebird_species_name)

normalize <- function(x) (x - min(x)) / (max(x) - min(x)) * 5 + 1  # Scale to range [1, 6]
species_urban_response$beak_width_norm <- normalize(species_urban_response$beak_width)
species_urban_response$hwi_norm <- normalize(species_urban_response$hwi)
species_urban_response$mass_norm <- normalize(species_urban_response$mass)

g = ggtree(tree_pruned, layout="fan") %<+% species_urban_response + abundance_proxy_scale

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

g + 
  geom_fruit(geom = geom_point, mapping = aes(size = beak_width_norm), fill = beak_width_colour, colour = 'black', shape = 21, stroke = 0.2, offset = 0, show.legend = FALSE) + 
  geom_fruit(geom = geom_point, mapping = aes(size = hwi_norm), fill = hwi_colour, colour = 'black', shape = 21, stroke = 0.2, offset = 0.13, show.legend = FALSE) + 
  geom_fruit(geom = geom_point, mapping = aes(size = mass_norm), fill = mass_colour, colour = 'black', shape = 21, stroke = 0.2, offset = 0.15, show.legend = FALSE) + 
  geom_tiplab(size=2, aes(label=ebird_species_name_clean, color = mean_relative_abundance_proxy), offset = 15) +
  labs(color = "Mean relative abundance proxy") + 
  theme(legend.position = "bottom") + ggplot2::xlim(0, 80)

ggsave(filename(FIGURES_OUTPUT_DIR, 'phylogeny.jpg'), width = 2000, height = 2100, units = 'px')

Format plot for nature

colors = c('Beak Width' = beak_width_colour, 'HWI' = hwi_colour, 'Mass' = mass_colour)
g_fig_base = ggtree(tree_pruned, layout="fan") %<+% species_urban_response

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

g_fig = g_fig_base + 
  geom_fruit(geom = geom_point, mapping = aes(size = beak_width_norm, fill = 'Beak Width'), colour = 'black', shape = 21, stroke = 0.2, offset = 0, show.legend = c(size = FALSE, fill = TRUE)) + 
  geom_fruit(geom = geom_point, mapping = aes(size = hwi_norm, fill = 'HWI'), colour = 'black', shape = 21, stroke = 0.2, offset = 0.13, show.legend = c(size = FALSE, fill = TRUE)) + 
  geom_fruit(geom = geom_point, mapping = aes(size = mass_norm, fill = 'Mass'), colour = 'black', shape = 21, stroke = 0.2, offset = 0.15, show.legend = c(size = FALSE, fill = TRUE)) + 
  geom_tiplab(size=2, aes(label=ebird_species_name_clean), offset = 30) +
  scale_fill_manual(name = "Relative trait size", values = colors) + 
  guides(fill = guide_legend(override.aes = list(size = 6))) + new_scale_fill()

heatmap_fig_ata = data.frame(species_urban_response[, "mean_relative_abundance_proxy"])
colnames(heatmap_fig_ata) = c('MRAP')
rownames(heatmap_fig_ata) = species_urban_response$ebird_species_name

gheatmap(g_fig, 
         heatmap_fig_ata, 
         offset = 14, 
         width = 0.2,
         colnames_angle = 0, 
         colnames_offset_x = 0.5, 
         font.size = 3) + 
  scale_fill_continuous(low = cyan, high = orange, guide = guide_colorbar(barwidth = 11)) + 
  labs(fill = "MRAP\n(Mean relative abundance proxy)") + 
  theme(legend.position = "bottom", legend.direction = "horizontal", legend.box = "vertical", legend.title.position = "top", legend.justification = "centre") + 
  ggplot2::xlim(0, 85)

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

ggsave(filename(FIGURES_OUTPUT_DIR, 'figure3.jpg'), width = 180, height = 220, units = 'mm', dpi = 450)

Plot traits against abundance proxy

cut_off_3 = 50
cut_off_2 = 35
cut_off_1 = 1

hwi_by_mass_urban_3 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_3) %>% slice(chull(hwi, mass))
hwi_by_mass_urban_2 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_2) %>% slice(chull(hwi, mass))
hwi_by_mass_urban_1 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_1) %>% slice(chull(hwi, mass))

beak_width_by_mass_urban_3 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_3) %>% slice(chull(beak_width, mass))
beak_width_by_mass_urban_2 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_2) %>% slice(chull(beak_width, mass))
beak_width_by_mass_urban_1 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_1) %>% slice(chull(beak_width, mass))

hwi_by_beak_width_urban_3 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_3) %>% slice(chull(hwi, beak_width))
hwi_by_beak_width_urban_2 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_2) %>% slice(chull(hwi, beak_width))
hwi_by_beak_width_urban_1 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_1) %>% slice(chull(hwi, beak_width))

species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_3) %>% nrow()
species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_2) %>% nrow()
species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_1) %>% nrow()

polygon_line_type = 'dashed'
polygon_linewidth = 0.4
polygon_alpha1 = 0.05
polygon_alpha2 = 0.1
polygon_alpha3 = 0.15

hwi_by_mass = ggplot(species_urban_response, aes(x = hwi, y = mass, colour = mean_relative_abundance_proxy)) + 
  geom_polygon(data = hwi_by_mass_urban_1, alpha = polygon_alpha1, color = trait_polygon_light, fill = trait_polygon_light, linewidth = polygon_linewidth, linetype = polygon_line_type) +
  geom_polygon(data = hwi_by_mass_urban_2, alpha = polygon_alpha2, color = trait_polygon_med, fill = trait_polygon_med, linewidth = polygon_linewidth, linetype = polygon_line_type) +
  geom_polygon(data = hwi_by_mass_urban_3, alpha = polygon_alpha3, color = trait_polygon_dark, fill = trait_polygon_dark, linewidth = polygon_linewidth, linetype = polygon_line_type) +
  geom_point() + 
  abundance_proxy_scale +
  labs(colour = "Mean relative abundance proxy", y = 'Mass', x = '') + 
  theme(legend.position = 'bottom', legend.title.position = 'top', legend.key.width = unit(12, 'mm'))

beak_width_by_mass = ggplot(species_urban_response, aes(x = beak_width, y = mass, colour = mean_relative_abundance_proxy)) + 
  geom_polygon(data = beak_width_by_mass_urban_1, alpha = polygon_alpha1, color = trait_polygon_light, fill = trait_polygon_light, linewidth = polygon_linewidth, linetype = polygon_line_type) +
  geom_polygon(data = beak_width_by_mass_urban_2, alpha = polygon_alpha2, color = trait_polygon_med, fill = trait_polygon_med, linewidth = polygon_linewidth, linetype = polygon_line_type) +
  geom_polygon(data = beak_width_by_mass_urban_3, alpha = polygon_alpha3, color = trait_polygon_dark, fill = trait_polygon_dark, linewidth = polygon_linewidth, linetype = polygon_line_type) +
  geom_point() + 
  abundance_proxy_scale + 
  labs(y = '', x = 'Beak width')

hwi_by_beak_width = ggplot(species_urban_response, aes(x = hwi, y = beak_width, colour = mean_relative_abundance_proxy)) + 
  geom_polygon(data = hwi_by_beak_width_urban_1, alpha = polygon_alpha1, color = trait_polygon_light, fill = trait_polygon_light, linewidth = polygon_linewidth, linetype = polygon_line_type) +
  geom_polygon(data = hwi_by_beak_width_urban_2, alpha = polygon_alpha2, color = trait_polygon_med, fill = trait_polygon_med, linewidth = polygon_linewidth, linetype = polygon_line_type) +
  geom_polygon(data = hwi_by_beak_width_urban_3, alpha = polygon_alpha3, color = trait_polygon_dark, fill = trait_polygon_dark, linewidth = polygon_linewidth, linetype = polygon_line_type) +
  geom_point() + 
  scale_fill_continuous(low = cyan, high = orange, guide = guide_colorbar(barwidth = 100)) + 
  labs(y = 'Beak width', x = 'HWI')

legend <- ggpubr::get_legend(
  # create some space to the left of the legend
  hwi_by_mass
)

plot_grid(
  nrow = 2, ncol = 2,
  hwi_by_mass + theme(legend.position="none"),
  beak_width_by_mass + theme(legend.position="none"),
  hwi_by_beak_width + theme(legend.position="none"),
  legend
)

ggsave(filename(FIGURES_OUTPUT_DIR, 'traits.jpg'), width = 2000, height = 1800, units = 'px')